Singular Behavior of Anderson Localized Wavefunctions for a Two-Site Model 
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We show analytically that the apparent non-analyticity discovered recently in the inverse par- 
ticipation ratio (IPR) of the eigenstates in Anderson's model of localization is also present in a 
simple two-site model, along with a concurrent non-analyticity in the density of states (DOS) at the 
same energy. We demonstrate its evolution from two sites to the thermodynamic limit by numerical 
methods. For the two site model, non-analyticity in higher derivatives of the DOS and IPR is also 
proven to exist for all bounded distributions of disorder. 
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I. INTRODUCTION 

In a recent paper [I], we reported numerical studies 
of the Anderson model of localization |2j in the insulat- 
ing regime. We studied, in particular, the ensemble- 
averaged inverse participation ratio (IPR), 1(E), which 
characterizes the extent of eigenstates at a given energy 
E. We found that 1(E) exhibits a sharp, apparently 
singular, behavior as a function of E, at a specific en- 
ergy E, which lies in the insulating phase, but is clearly 
distinct from the true band edge. This sudden change 
in behavior was interpreted in terms of a sharp transi- 
tion from a regime comprising for the most part of typi- 
cal, Anderson-localized states, to one comprising almost 
exclusively states involving resonance between two (or 
more) sites. Such states involving resonance among a 
large number of sites are known in the literature as Lif- 
shitz states [3J, [J- 

This aforementioned behavior was seen in one, two and 
three dimensions for bounded distributions of on-site en- 
ergies (e.g. the uniform distribution of width w with a 
constant probability density between the limits ±w/2.), 
but not for unbounded distributions (e.g. Gaussian). 
The apparent singular behavior could be observed most 
clearly for disorder exceeding a certain, moderate value 
(disorder width w — 3.8 for the uniform distribution in 
d = 1), and appeared to persist for extremely large dis- 
order, where most states at the energies of interest are 
localized on a few sites. This provides motivation for 
studying various quantities analytically for the Anderson 
model defined on a small, finite number of sites. The 
simplest of these exhibiting the difference between typi- 
cally localized and resonant states is the two-site model, 
which, as we show below, can be solved analytically. 

The solution for the two-site model shows that 
for all bounded distributions defined in the interval 
(—w/2,w/2), both the ensemble averaged IPR, 1(E), 
and the density of states p(E), exhibit singular behav- 
ior at the energies E = ±E 2 whose value is related to 
the diso rder width w by the simple algebraic expression 
E% = -y/u> 2 /4 + 1. The exact nature of the singularity 



depends on the form of the distribution characterizing 
the diagonal disorder of the Anderson model. In partic- 
ular for the uniform distribution of Anderson's original 
paper[5], with a step discontinuity at —w/2 and +w/2, 
both quantities show a slope discontinuity. This is rather 
similar to what is seen numerically for the thermody- 
namic system PQ, though the energy at which the appar- 
ent singular behavior is seen for the thermodynamic sys- 
tem differs from that of the two-site model. Nevertheless, 
the reason for the discontinuity in slope is found to be 
the same in the two-site model, namely the loss of cer- 
tain kind of states beyond a critical energy [T], so it is not 
unreasonable to conjecture that the two phenomena are 
related. 

For other bounded distributions, with power-law 
thresholds, the singularity for the two-site model appears 
in a higher order derivative, while for integrable, inverse 
power law thresholds (where the distribution diverges at 
the edges), the derivatives of p(E) and 1(E) diverge at 
E = ±E 2 . 

The plan of the paper is as follows. In Section II, we 
define the n-site Anderson model on a ring (i.e., with 
periodic boundary conditions). Since all states of the 
Anderson model are known to be localized in d = 1, this 
set of models are representative of the localized phase, 
and also of the true thermodynamic model for sufficiently 
large disorder. We then specialize to the case of two sites 
and present complete results for the ensemble averaged 
density of states and the inverse participation ratio for 
the case of Anderson's uniform distribution. The details 
of the calculation, provided in Appendix A, make it clear 
that a bounded distribution is necessary to obtain singu- 
lar points for the two-site model. We then consider other 
bounded distributions with power law singularities, and 
demonstrate the singular behavior of appropriate deriva- 
tives of p(E) and 1(E). Details of the derivation are 
provided in Appendix B. In Section III, we present our 
results for p(E) and 1(E) for the n-site Anderson model 
with a uniform distribution obtained numerically for two 
values of w, representing moderate and large disorder, 
and show how these quantities evolve with n. Finally, in 
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The IPR for a wavefunction "J = J2i a i\i > is defined 



FIG. 1: Some normalized distributions of the on-site energy 
considered in this paper. On the left are three bounded dis- 
tributions: (i) Uniform, P(e) = 1/w, (ii) & (iii) have power 
law edges, P(e) = Jjffffiff, (1 - 4e 2 /w 2 ) x for A = 1/2 and 
A = —1/2 respectively. On the right is the unbounded Gaus- 
sian distribution. 



Section IV, we summarize our results. 



II. ANDERSON MODEL AND SYMMETRIC, 
BOUNDED PROBABILITY DISTRIBUTIONS 

The Anderson model Hamiltonian is given by: 

H = ^2ei\i><i\ + VY;\i><j\ (1) 

i <i,j> 

where \i > are states localized on sites i of a simple hyper- 
cubic lattice in d-dimensions, and i, j are nearest neigh- 
bours. The onsite energies q are independent random 
variables taken from a specified distribution P(e). By 
definition, J™ P(e)de = 1. We take P(e) to be non-zero 
only between ±w/2. Within these limits, the distribution 
can have any form as long as it is continuous and sym- 
metric in energy. Examples are the uniform, semicircu- 
lar and inverse-power law distributions (Fig[l](i) - (iii)). 
Later in this section, we will briefly discuss what happens 
when P(e) is not bounded (e.g. Gaussian distribution 
(FigfTl (iv))). In the thermodynamic limit, for symmet- 
ric P[e) and d-dimensional hypercubic lattices which are 
bipartite (including the nearest neighbour chain that we 
consider here), properties are symmetric around E = 0. 
The symmetry in energy also holds for finite size systems 
that have an even number of sites. 

The hopping V is assumed to have the same non-zero 
value for all nearest neighbours, and zero for non-nearest 
neighbours. In this paper, we use V = — 1 to define our 
unit of energy. In one-dimension, which is what we con- 
sider in this paper, all eigenstates in the Anderson model 
are known to be localized [5] for non-zero w. The local- 
ization length decreases as the disorder, parametrized by 
the range of possible values of e< (e.g. the width of its 
probability distribution w), increases. 
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It is inversely proportional to the width (support) of the 
wavefunction. For localized eigenstates, the IPR goes 
to a constant as the system size tends to infinity. As 
detailed in our previous paper, in the Anderson model 
which uses a uniform distribution, the IPR has a local 
(non-zero) minimum at the centre of the band, rises to 
a maximum and turns around sharply in an apparently 
non-analytic manner and goes to zero at the band edges 
E = ±{w/2 + 2). 

We now consider the Anderson model for finite lattices. 
For two sites (1,2), the Hamiltonian of Eq. [Tjreduces to 

H = e 1 \l >< 1| + e 2 |2 >< 2| - |1 >< 2| - |2 >< 1| (3) 

i.e. a 2 x 2 matrix, which one can easily solve for arbitrary 
ei and e 2 . The eigenvalues are 



where x = £l J^ 2 and y 



£l-£2 
V2 



(4) 



The corresponding 



eigenstates can be written as a\l > +/3\2 > where 




The IPR for both wavefunctions is the same: 
I + =l- 1 



1 + 2/V2 



(5) 



(6) 



which is further independent of x. 

The random variables e± and e 2 have a probability dis- 
tribution P(e). The ensemble averaged DOS (p) and IPR 
(J) are given by: 

p{E) = - / de 1 de 2 P(e 1 )P(e 2 ) x 

^ J -oo J-oo 

[6{E - Ex) + 6{E - E 2 )] (7) 



1(E) 



1 1 



/oo poo 
/ de 1 de 2 P(e 1 )P(e 2 ) x 
-oo J —oo 



p{E) 2 

[6(E - Ex) + 6{E - E 2 )] x 
1 



1 - 



2(1 + f) 



(8) 



Transforming the variables to x and y makes the integrals 
easier to evaluate. The Jacobian is unity, giving dt\dt 2 = 
dxdy. 

For w — 0, the DOS will be two delta functions at 
E = ±1. For < w < 2, the DOS goes to zero for 
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V2 x = £ i + £ 2 V2 

V2 

FIG. 2: Allowed values of x and y lie inside the square. Color- 
bar shows variation of IPR over square. The red curve shows 
the left branch of the hyperbola when it is tangent to sides 1 
and 4 of the square. 



w/2 - 1 < E < -w/2 + 1 (see Appendix A). Since the 
two-site model is expected to be a good approximation 
for the infinite model only at large w, we take w > 2 from 
now on. 

Under these conditions, the integrands in Eq. [7] and 
Eq. [8] are non-zero only over the square shown in Fig.l. 
As shown in Appendix A, the arguments of the delta 
functions in the integrals are just two branches of a hyper- 
bola symmetric about the x-axis, (x/\/2—E) 2 — y 2 /2 = 1. 
At E = 0, it is also symmetric about the y-axis. As E 
increases (decreases), the hyperbola moves right (left) as 
shown in Fig [7] of Appendix A. Any ensemble averaged 
property of the wavefunctions can be calculated as a func- 
tion of energy by doing the relevant line integral along 
the length of the hyperbola lying inside the square. 

The simplest bounded probability distribution is one 
in which the P(e) is a constant between —w/2 and w/2. 
This is also the distribution considered in Anderson's 
original paper[2]. The complete expressions for the DOS 
and IPR when the probability distribution is uniform are 
given in Table [I] and the derivation is provided in Ap- 
pendix A. In Fig [3] and Fig [IJ we plot (as solid lines) 
the analytical DOS and IPR for 2 sites as a function of 
energy for w = 4.5 and w = 6 respectively which agree 
with our numerical results shown as rightward-pointing 
triangles. Going from the centre to the edge of the band, 
the slope of the DOS and IPR discontinuously changes 
sign at E = ±E%. This can be understood as a loss of 
typical "Anderson" type states as follows. 

The IPR is an average over all the states lying on the 
hyperbola. Each point in the square will contribute to 
the average at two different energies, once when it lies 
on the left branch of the hyperbola and once when it 
lies on the right. The states with the highest IPR lie on 
the upper and lower vertices of the square. These states 
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TABLE I: Analytical expressions for DOS and IPR over the 
whole band for the uniform distribution when w > 2. Here, 
a = E + w/2, b = E - w/2 and a = ^fi, /3 = 



occur when \y\ is maximum, i.e. t\ and e 2 are far from 
each other. Near the center of the square, ei is close to 
e 2 — wavefunctions arising from these configurations can 
be called 'resonant' states. 

At E = E2, the left branch of the hyperbola becomes 
tangent to the sides 1 and 4 of the square as shown in Fig. 
1. In this position, the maximum length of the hyperbola 
passes through the high IPR region. The average IPR 
will thus be maximum in this region. However, at any 
energy close to this, the hyperbola shifts and the high- 
IPR states start being lost leading to a decrease in the 
average IPR. The rate of decrease is different when the 
shift is left or right, leading to a discontinuity in the slope 
dI }^ of the average IPR. When the hyperbola moves left 
(E decreases), it still passes through a large number of 
high IPR states. When it moves right (E increases), it 
loses the high-IPR states much faster, and passes through 
states for which t\ + e 2 is large and t\ — £2 is small, im- 
plying that both sites have energies close to w/2. These 
type of resonant states which occur near the edge of the 
band and arise from clusters of sites which have energies 
close to one of the band edges, are called Lifshitz states. 
They are separated here from the rest of the band by the 
non-analytic point in the IPR. 

The softer discontinuity in the DOS is only due to a 
decrease in the total number of states lying on the hyper- 
bola, irrespective of whether they are regular Anderson 
localized or resonant states. 

For a non-constant (but still bounded) probability dis- 
tribution, the discontinuity may be present in higher 
derivatives of the DOS and IPR. As shown in Appendix 
B, the expression for (d n p/dE n ) or (d n I/dE n ) will have 
terms of the form (d mi P(e 1 )/de™ 1 )(d m2 P(e 2 )/rfe 2 ' 12 ) 
with (mi + ^2) ranging from to (n — 1). If any of 
these terms are non-zero at the edge, there will be a dis- 
continuity in that derivative at E = ±E%. If any of them 
diverge at the edge, then the corresponding (d n p/dE n ) 
and (d n I/dE n ) will also diverge at E = ±E 2 . Thus, if a 
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probability distribution P(e) — > (1— Ae 2 /w ) as e — !• w/2 
with — 1 < A < 0, then all the derivatives of DOS and 
IPR will be infinite at E = ±E 2 - For A positive, the 
divergence will be in a higher derivative. 

From the above discussion, it is easy to see that a non- 
bounded probability distribution like the Gaussian, for 
which the boundaries of the square lie at infinity, will 
not exhibit any discontinuities in the 2-site model. 



III. NUMERICAL RESULTS 



In Fig [3] and Fig [4j we plot with varying symbols nu- 
merical results for n-site chains for typical values of dis- 
order, w = 4.5 and 6. These typical cases show how the 
DOS and IPR converge to their thermodynamic limits 
relatively quickly as n is increased, as expected since we 
are in the localized regime. Larger values of w (larger 
disorder) require smaller size systems to reach the ther- 
modynamic limit within a given precision as the localiza- 
tion length decreases with increasing disorder. Thus, for 
W = 4.5, a lattice size of 100 is required to obtain results 
significantly closer to the thermodynamic limit than the 
size of the symbols, whereas for w = 6, only 20 sites 
are required. The discontinuity in the derivative of the 
IPR survives and seems to grow sharper in the thermo- 
dynamic limit, whereas that for the DOS seems to flatten 
out a bit. There also seems to be an additional feature in 
the IPR beyond the maximum. This may be the start of 
a transition into another regime of resonant states, e.g. 
from 2 sites to 3 sites. However, this is hard to con- 
firm numerically. Interestingly, as would be expected if 
the above conjecture were true, it is not present in the 
two-site model. 

To measure the evolution of the non-analytic behav- 
ior with size, the curves on either side of the roughly 
estimated maxima of the DOS and IPR were fit with 
quadratic functions up to the points where additional 
kinks occur in the curve. The position of the maxima, 
Ep.i and Ej_\ were then determined from the intersec- 
tion of these fits. The difference of the slopes si e f t and 
Sright of the fits near the maxima were used to quantify 
the discontinuity. (As = s; e /t — s r i g ht) is shown to reach 
a constant value with increasing size in Fig [5] suggesting 
that the singularity survives in the thermodynamic limit. 

Fig [6] compares the analytically determined point of 
discontinuity (E = y // w 2 /4 + l)in the 2-site case to the 
numerically determined one in the infinite size limit. For 
large disorder, the 2-site problem seems to be a good ap- 
proximation to the infinite size limit, whereas for smaller 
w the deviation is larger. This is to be expected since 
increasing disorder leads to smaller localization lengths. 




FIG. 3: Evolution of DOS and IPR with increasing system 
size for w = 4.5. Solid line shows analytical result for n—2. 



IV. CONCLUSIONS 

In conclusion, we have shown by direct calculation how 
the non-analyticity in the DOS and IPR arises for the 2- 
site model for bounded disorder. Further, we have shown 
the evolution with number of sites for the case of uniform 
disorder and found an apparent singularity, at least for 
IPR, surviving in the thermodynamic limit. Our numer- 
ical studies on different lattice sizes suggest that there is 
an exact energy at which high IPR states are lost and 
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FIG. 4: Evolution of DOS and IPR with increasing system 
size for w — 6. Solid line shows analytical result for n=2. 



that the position as well as the magnitude of the singu- 
larity approach a limit as n —¥ oo. It is at least clear 
that the singularity arises due to the bounded nature 
of the disorder distribution, which leads to loss of high- 
IPR wavefunctions at a particular energy in the band like 
in the 2-site model. The non-analyticity appears to be 
present for all lattice sizes and we surmise that it survives 
in the thermodynamic limit. The question of whether it 
is possible to show that the non-analyticity persists in 
the thermodynamic limit by analytical means remains 
unanswered at present. 
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FIG. 5: As = (sieft — s r ight) as a function of system size for 
w = 6 for DOS (triangles) and IPR (circles). 
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FIG. 6: Comparison of value of energy at which singularity 
occurs for 2 sites (blue solid line) and the large size limit 
(red circles). The red dotted line is a fit to the circles. The 
leading term w/2 is subtracted from the energy to amplify 
the difference. 



We also remark that studying a 3-site model may allow 
us to determine if the shoulders seen in the IPR curves 
in Fig. [3] and Fig. [4] are due to transition from 2 to 3 
resonant sites. 

After this work was completed, we received news of 
parallel work by Ujfalusi & Varga [5] which addresses 
some of the issues discussed in this paper. 

This work was supported by DOE grant de- 
sc20002140. 



V. APPENDIX A 

The integrands in Eq. [7] and Eq. [8] contain delta func- 
tions in a 2 dimensional space. Using the identity for n- 



G 



dimensions: 



8(g[r))f{r)df 



da(r) 
g(f)=o |Vff| 



(9) 



where the integral on the right is over the hypersurface 
g(f) = 0, in one dimension lower than the one on the left. 
In our case g(f) = E— A^(x± \Jy 2 + 2) i.e. the equation 

of a hyperbola for fixed g(f), here g(r) — 0. When E = 0, 
the hyperbola is symmetric about the y-axis, and when E 
is non-zero, the hyperbola is displaced along the s-axis. 
It is always symmetric about the x-axis. 

Because of the delta function and the bounds on x and 
y, we effectively have to calculate a line integral along 
the portion of the hyperbola which lies inside the square. 
The hyperbola can be parametrized by x — \/2{sec(t) + 

E), y = v^tan(t). Thus, we have |Vp| = = 



2 2^(0 1 and da = \/ 2 sec 2 (t)t&n 2 (t) + 2 sec 4 (t)dt = 



•\/2sec(i)-\/tan 2 (t) + sec 2 (/). The intersection points of 
the sides of the square (numbered as in Fig.l) with the 

(f+E) 2 -l (%-E) 2 -l 



V2(f-E) ' 



yi 



hyperbola are: y 3 = *^ {m _ +E) , 2/4 

-3/4 and y 2 = -J/3- 

The five possible situations for w > 2 are illustrated in 
Fig. [7] The size of the square depends only on w. The 
position of the hyperbola depends only on the energy. For 
w < 2, the density of states goes to zero in the middle 
case (Fig. [7]hdi)) since the square is small enough to 
lie between the two branches of the hyperbola without 
intersecting either of them. 

In cases (i) and (ii), only the right branch of the hy- 
perbola (E-) falls inside the square, for (iv) and (v), we 
need to account for only the left branch (E + ), and in case 
(hi), both branches are required. For positive energies, 
the point which gives rise to the change in sign of the 
slopes of average DOS and IPR is at the transition from 
case (iv) to case (v) (There will also be a smaller discon- 
tinuous change in the slopes when transitioning from (hi) 
to (iv)). 

The integrand is independent of E and the dependence 
on E is provided by the limits of the integration. The 
average IPR and DOS for the uniform distribution are 



p(E) 



yo(E) 



Vo(E) 



2 sec 2 (i) eft 



Similarly, 
1(E) = 



yo(E) 



P(E) J-y o{E ) 



2sec 2 (i) 1 



2sec 2 (t) 



(10) 



dt (11) 



with limits of the above integrals given by the intersection 
points and using one or both branches of the hyperbola 
as appropriate. The final expressions are given in Table 

m 

The transitions between case (i) to case (ii) and case 
(iv) to case (v) occur at E = ±2?2- At these points, 




(iii) 



(iv) 



(x±^/y 2 + 2) moves to the right as 



FIG. 7: Hyperbola E 
E increases. Figures (i)-(v) depict sequentially the five cases 
in Table [I] Allowed values of x and y lie inside the square 



2/3 = 2/4 = 



\/2- 



From Eq. 2, x — at his point. This 

corresponds to the situation t\ = w/2 and e 2 = —w/2 or 
vice versa, that is, the two sites lie at opposite ends of 
the probability distribution. 



VI. APPENDIX B 

The DOS (and IPR) contains contributions p\ and p 2 
from the 2 branches of the hyperbola. In this section, the 
analysis is carried out only for p\. The other term can 
be calculated in an analogous manner. 

For any symmetric, bounded and continuous probabil- 
ity distribution P(e), the DOS due to the left branch of 
the hyperbola can be expressed in x and y coordinates 
as: 




p [ *( E >v) + v\ , 



V2 



(12) 
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In the integrand x is now a function of E and y: x = 

V2E + P (^f ) = P(ei) and P (^) = 

P(e2). The limits yo(E) depend on the width of the prob- 
ability distribution. 

Differentiating Eq{l2]with respect to E, 



dpi(E) 
dE 



V2 



^ dy d(P(e 1 )P(e 2 )) 



dE 



[P(ex)P(e 2 )] l 



Consider the left and right limits of 



dyo 
1 dE 



dpi(E) 
dE 



(13) 



at E = 



^- + 1 when the probability distribution is bounded 
between —w/2 and w/2. yo = 2/3 (intersection of the 
hyperbola with sides 2 and 3 of the square) for the left 
limit and yo = 2/4 (intersection with sides 1 and 4) for 
the right limit. The first term (the integral) remains the 

same for both since at the critical energy E = — + 1, 

j/3 = j/4 = w/\j2. Now consider the second term. 
% jk P( ei ) = P( w /2) and P(e 2 ) = P(-w/2), 

i.e. equal to their value at the edges. For the uniform 
distribution, these are non-zero, and therefore the first 
derivative of DOS is discontinuous. 

Similarly, we can calculate the second derivative: 



and so on. Note that 



dP(e) de 
de dE' 



The integration term on the right is always same 
for both ^limits. fc The other fc terms are of the form 
d defi*^ d de™* 2) ^d~M- Since 1W is alwa y s discontinu- 
ous at E — Ei for any k 7 a non-zero or infinite value of 
d cL"*^ 1 " 1 d de™^ i s required for a discontinuity in p. For 
example, if P(e) has the form (1 — 2x/w) x at the edge, 
then 2A — mi — rri2 > for analyticity. Therefore, the 
derivative on the left side of the equation will be discon- 
tinuous at E = E-i if it is of order (n + 1) or greater where 
n is the smallest integer greater than or equal to 2A. 



d 2 Pl (E) 
dE 2 



V2 



Vo dy d 2 (P(e 1 )P(e 2 ))\ 



o dE* 
d(P(e 1 )P(e 2 )) y=yo 
dE 

d(P( £l )P(6 2 )) 

dE 



(P(e 1 )P(e 2 )) y=ya 



y=vo 
d 2 y 
dE 2 



dyo 
dE 

dyo 
dE 



(14) 



For a non-bounded distribution, yo will be infinity, in- 
dependent of E and therefore the DOS is always smooth. 
In the 2-site problem, since the IPR and its derivatives 
are only a function of y and independent of x, they will 
have discontinuities at the same orders and energies as 
the DOS. 
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